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Abstract. We calculate the ionisation fraction in protostellar disk models using two different gas-phase chemical 
networks, and examine the effect of turbulent mixing by modelling the diffusion of chemical species vertically 
through the disk. The aim is to determine in which regions of the disk gas can couple to a magnetic field and 
sustain MHD turbulence. The disk models are conventional Q-disks, and the primary source of ionisation is X-ray 
1 irradiation from the central star. We assume that the vertical mixing arises because of turbulent diffusion, and 

accordingly equate the vertical diffusion coefficient, T>, with the kinematic viscosity, v. 

We find that the effect of diffusion depends crucially on the elemental abundance of heavy metals (magnesium) 
included in the chemical model. In the absence of heavy metals, diffusion has essentially no effect on the ionisation 
structure of the disks, as the recombination time scale is much shorter than the turbulent diffusion time scale. 
When metals are included with an elemental abundance above a threshold value, the diffusion can dramatically 
reduce the size of the magnetically decoupled region ("dead zone"), or even remove it altogther. This arises when 
recombination is dominated by metal ions, and the recombination time exceeds the vertical diffusion time scale. 
For a complex chemistry the elemental abundance of magnesium required to remove the dead zone is iM g = 10 -10 - 
\f~^ ' 10~ 8 . We also find that diffusion can modify the reaction pathways, giving rise to dominant species when diffusion 

is switched on that are minor species when diffusion is absent. This suggests that there may be chemical signatures 
of diffusive mixing that could be used to indirectly detect turbulent activity in protoplanetary disks. 
We find examples of models in which the dead zone in the outer disk region is rendered deeper when diffusion 
is switched on. This is caused by turbulent mixing diluting the electron fraction in regions where the ionisation 
degree is marginally above that required for good coupling. 

Overall these results suggest that global MHD turbulence in protoplanetary disks may be self-sustaining under 
favourable circumstances, as turbulent mixing can help maintain the ionisation fraction above that necessary to 
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ensure good coupling between the gas and magnetic field. 

Key words, accretion, accretion disks - MHD - planetary systems: protoplanetary disks - stars: pre-main sequence 

1. Introduction cosity. At the present time the only source of turbulence 

in accretion disks that is known to work is the magnetoro- 

Observational studies of young stars in star forming re- tational instability (MRI), which gives rise to MHD tur- 

gions indicate that protostellar disks are a common occur- bulence (Balbus & Hawley 1991; Hawley & Balbus 1991; 

rence (e.g. Beckwith & Sargent 1996; Prosser et al. 1994). Hawley, Gammie & Balbus 1996). 

Many disks show signatures of active accretion, with the As has been well documented in the literature, 

probability of accreting disks being present, and the app- there are continuing questions about the applicability 

parent gas accretion rate, scaling inversely with the age of of the MRI to protostellar disks because of their high 

the stellar system in which the young stars are embedded, densities and low temperatures that lead to low lev- 

The canonical value for the gas accretion rate, however, els of ionisation (e.g. Blaes & Balbus 1994; Gammie 

is often quoted as being M ~ 10~ 8 M Q yr _1 (e.g. Sicilia- 1996). Magnetohydrodynamic simulations of disks includ- 

Aguilar et al. 2004). ing ohmic resistivity (Fleming, Stone & Hawley 2000) 

There are a number of potential mechansims that may show that for magnetic Reynolds numbers Re m smaller 

lead to angular momentum transport in protostellar disks, than a critical value -Re™*, turbulence cannot be sustained 

giving rise to the observed mass accretion. Angular mo- and the disks return to a near-laminar state. Typically the 

mentum transport that occurs globally throughout the gas-phase electron fraction should be x[e~] ~ 10 -12 for 

disk, producing accretion at the observed rates, probably disks to be able to sustain MHD turbulence, 
requires turbulence to act as a source of anomalous vis- 
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A number of studies of the ionisation fraction in pro- 
tostcllar disks have appeared in the literature. Gammie 

(1996) suggested that disks may have magnetically "ac- 
tive zones" sustained by thermal or cosmic ray ionisation, 
adjoining regions that were magnetically inactive - "dead 
zones". More recent studies have examined this issue in 
greater depth. Sano et al. (2000) used a more complex 
chemical model that included dust grains. Glassgold et al. 

(1997) and Igea et al. (1999) examined the role of X-rays 
as a source of ionisation in protoplanetary disks. Fromang 
et al. (2002) considered the role of heavy metals in de- 
termining the ionisation fraction because of the poten- 
tial importance of charge-transfer reactions. Matsumura 
& Pudritz (2003) examined the ionisation fraction in the 
externally heated, passive disk model proposed by Chiang 
& Goldreich (1997) using the Sano et al. (2000) chemical 
reaction network. Semenov et al. (2004) recently studied 
disk chemistry using a complex reaction set drawn from 
the UMIST database. 

In a recent paper (Ilgner & Nelson 2005 - hereafter 
"paper 1"), we examined and compared the predictions 
made by a number of chemical reaction networks about 
the ionisation fraction in standard a-disk models. This 
study included an examination of the reaction scheme pro- 
posed by Oppenheimer & Dalgarno (1974), and more com- 
plex schemes drawn from the UMIST database, in addi- 
tion to a number of gas-grain chemical networks. In this 
paper we extend this initial study to examine the role of 
turbulent mixing in determining the ionisation fraction 
in protoplanetary disk models using gas-phase chemistry. 
We consider the simple Oppenheimer & Dalgarno (1974) 
model, and a more sophisticated chemical reaction net- 
work based on the UMIST database. We allow the ver- 
tical diffusion of each chemical species to occur, equat- 
ing the diffusion coefficient T> to the kinematic viscosity 
v that drives disk accretion. We find that the inclusion of 
diffusion can have a very significant effect on the ionisa- 
tion structure of protoplanetary disks, in particular when 
a small abundance of heavy metals (magnesium) is in- 
troduced into the reaction networks. In some cases the 
disk can modify the ionisation fraction sufficiently that 
the dead zone disappears entirely. 

This paper is organised as follows. In section [21 we de- 
scribe the chemical models, our implementation of the dif- 
fusion equation for chemical species, and the numerical 
method used to solve the reaction-diffusion equations. In 
section we discuss the various models that we compute 
and present their results. In section 01 we summarise our 
main findings and discuss their potential consequences for 
protoplanetary disks. 

2. Reaction-diffusion model 

We consider a system of s chemical species whose local 
abundances in the protoplanetary disk evolve due to chem- 
ical reactions and vertical diffusive transport, arising be- 
cause of concentration gradients and driven by turbulence. 
In the discussion below we define r to be the number of 



chemical reactions, and n to be the total molar density 
(expressed in units of mols cm -3 ). The global elemental 
composition of the system is conserved by applying ap- 
propriate boundary conditions in our reaction-diffusion 
model. 

2.1. Disk models 

The underlying disk models considered are standard a- 
disks. Details are given in paper 1 and references therein. 
To recap: the disks are assumed to orbit a young solar 
mass star and undergo viscous evolution. We use the a 
prescription for the viscous stress, such that the kinematic 
viscosity v — aCg/Ct, where c s is the sound speed and £1 is 
the local Keplerian angular velocity. Heating of the disk is 
provided by viscous dissipation alone, and cooling by ra- 
diation transport in the vertical direction. The disk struc- 
ture is obtained by solving for hydrostatic and thermal 
equilibrium. We employ 30 zones in the vertical direction 
(from midplane to disk surface) and 100 zones in radius 
between 0.1 < R < 10 AU when computing the chemi- 
cal evolution. The underlying disk model was computed 
using 100 cells in the vertical direction, and values were 
interpolated onto the 30 grid points used in computing the 
chemistry. The disk models are completely specified by the 
mass accretion rate, M and the value of a. We consider 
two models: one with Al = 10~ 7 M Q yr _1 and a — 10 -2 , 
the other with M = 10~ 8 M Q yr -1 and a = 5 x 10~ 3 . The 
mass contained in these models is 0.0087 M Q and 0.0049 
M Q respectively. Contour plots showing the distribution 
of the kinematic viscosity, v, are shown in figures ^ and FJJ 

2.2. Kinetic model 

We have applied two different kinetic models for the gas- 
phase chemistry. The first is a simple five component 
model introduced originally by Oppenheimer & Dalgarno 
(1974). The second is a much more sophisticated reac- 
tion scheme based on the UMIST database (Le Teuff et 
al. 2000). We assume that ionisation of the disk mate- 
rial arises because of incident X-rays that originate in the 
corona of the central T Tauri star. We neglect contribu- 
tions from Galactic cosmic rays as they are not expected 
to penetrate into inner disk regions we consider due to 
the stellar wind. The details of these models have been 
described in paper 1, where they were given the labels 
modell and model3. We maintain this labelling scheme in 
this paper so as to establish continuity with our previous 
work. 

The kinetic equations obey the law of atomic balance: 

s 

J>y=a, Cj=1.2,...,r). (1) 

i=i 

Here i/y divided by the molecular mass M i of the i th com- 
ponent (defined in units of grams/mol) is proportional to 
the stoichiometric coefficient, while the i th species is in- 
volved in the j th chemical reaction. The molar density, n i , 
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Fig. 1. The kinematic viscosity v\ the contour lines refer 
to values 1 x 10 15 ,2 x 10 15 ,3 x 10 15 cm 2 s -1 . The disk 
parameters are a — 1CP 2 and M = 10 -7 MQyr" 1 . 



Fig. 2. The kinematic viscosity v\ the contour lines refer 
to values 3 x 10 14 ,4 x 10 14 ,5 x 10 14 cmV 1 . The disk 
parameters are q = 5x 1CP 3 and M = 10~ 8 M^yr -1 . 



of the i component is related to the mass density, p,, of 
the i th component through 

Pi = M ini . (2) 

The rate of change of the molar density of the i th compo- 
nent within a given volume due to chemical reactions and 
flow of species i into that volume is 

^ = -V.(n i v i ) + j2»i j Jr (3) 

Jj denotes the chemical reaction rate associated with 
the 7 th chemical reaction, while u i4 Jj is the forma- 
tion/destruction rate of the i th component due to the j th 
chemical reaction. Vj is the velocity of the i th component. 

2.3. Diffusion 

The first term on the right hand side of equation |[5J rep- 
resents the rate of flow of species i into a given volume. 
In this paper we assume that the underlying protoplane- 
tary disk is turbulent, and that turbulent diffusion acts to 
transport chemical species from one region of the disk to 
another. As we are primarily interested in how diffusion 
affects the ionisation fraction of the disk material, and 
we expect the largest gradients in the electron fraction to 
be in the vertical (z) direction, we consider only vertical 
diffusion in the paper. Equation becomes 

-ot = o- z { nV d- z x >) + ^ J » (* = 1 "--< s ) ( 4 ) 

\ / j=l 

where T> is the turbulent diffusion coefficient, and x% = 
rii/n is the fractional abundance of species i. We present 
a fuller discussion of the derivation of equation (J2J in ap- 
pendix. A similar set of equations has been used by nu- 
merous authors who have examined the effect of turbulent 



diffusion on the chemical state of molecular clouds (Xie, 
Allen & Langer 1995; Rawlings & Hartquist 1997; Willacy, 
Langer & Allen 2002) and on the role of radial mixing 
in protoplanetary disks (Wehrstedt & Gail 2002, 2003). 
In this work we adopt the approximation that D = i>, 
where v is the kinematic viscosity that drives the radial 
diffusion of mass through the protostellar accretion disk. 
Although there is still some debate about the precise rela- 
tionship between the vertical diffusion coefficient and the 
effective kinematic viscosity that arises because of MHD 
turbulence in astrophysical disks (see Carbillado, Stone & 
Pringle 2005; Johansen & Klahr 2005), it is not expected 
that T>/v differs greatly from unity. As already discussed, 
the basic kinetic models that we consider in this paper are 
models model 1 and model3 described in paper 1. When 
these kinetic models are used in conjunction with turbu- 
lent mixing we use the labels modellD and model3D, re- 
spectively. When diffusion is switched off we often use the 
term "pure kinetic models" when describing the results in 
later sections. 

2.4. Numerical method 

The reaction-diffusion model is governed by the set of s 
coupled parabolic partial differential equations 10} . These 
equations can be interpreted as a linear superposition of 
two operators, one describing the mass transfer of species 
i due to diffusion, and the other describing changes due 
to chemical reactions. There are number of possible ap- 
proaches to solving these equations. We employed the 
method of lines, which is a technique used, e.g., in atmo- 
spheric physics that solves the full system of equations si- 
multaneously (e.g. Chang et al. 1974), rather than a more 
usual operator splitting approach. Hence, we transformed 
the system of s PDEs into a system of s x nz ODEs, 
where nz is the number of grid points in the z direction. 
Adopting a uniform mesh along this coordinate direction, 
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we replaced the spatial derivative in equation by a fi- 
nite differencing scheme. 

The boundary conditions are taken to be symmetric at 
z = (disk midplane) and no flux at z = H (photospheric 
disk height) for all times t. Applying the Crank-Nicholsen 
finite differencing scheme and the chosen boundary condi- 
tions satisfy the conservation of elements and charges at 
any cylindrical radius R. Using nz = 30, we found that 
the elements and charges are conserved at all cylindrical 
radii and for all times t. At t < 10 5 yrs, for example, the 
change in both the elemental abundances and the total 
charge was below 10 -8 per cent compared with the corre- 
sponding values at t = yr. 

We have allowed all species, including molecular hy- 
drogen and helium, to be mixed by diffusion. We thus 
checked a posteriori if the assumed constancy of the mean 
molecular weight /i = 2.33 is satisfied. We find that it is, 
as expected from the fact that the value of \x is dominated 
by hydrogen and helium, and these species do not develop 
significant concentration gradients, thus precluding signif- 
icant diffusion. 

We also allowed free electrons to be mixed by diffu- 
sion. Differences in the concentration gradients of the ions 
and electrons can in principle lead to different diffusion 
velocities which may result in a non-zero diffusion flux of 
charges per zone. We found, however, that this did not 
occur. 

2.5. Initial conditions 

The reaction-diffusion calculations are initiated with all 
chemical elements being in neutral atomic form, apart 
from hydrogen which is assumed to be in molecular form. 
The species were distributed in the disk uniformly. Details 
are given in paper 1. Chemical changes lead to concentra- 
tion gradients that then initiate the action of diffusion. 



3. Results 

We have evolved the disk chemistry using the kinetic mod- 
els described in section l2~2l including mixing due to turbu- 
lent diffusion. The primary aim of this work is to compare 
and understand the differences in the distribution of free 
electrons that arise in the different models, and in partic- 
ular to understand the role that turbulent mixing plays 
in determining the free electron fraction and distribution. 
We wish to determine which parts of the disk are suffi- 
ciently ionised for the gas to be well coupled to the mag- 
netic field, and thus able to maintain MHD turbulence, 
and which regions are too neutral for such turbulence to 
be maintained. We refer to those regions as being "active" 
and "dead" zones respectively, with the region bordering 
the two being the "transition" zone. The important dis- 
criminant that determines whether the disk is active or 
dead is the magnetic Reynolds number, Re m , defined by 



where H is the disk semi-thickness, c s is the sound speed, 
and /i m is the magnetic diffusivity (not to be confused 
with the mean molecular weight). Numerical simulations 
(e.g. Fleming, Stone & Hawley 2000) indicate that a crit- 
ical value of the magnetic Reynolds number, i?e mt , ex- 
ists such that non linear MHD turbulence cannot be sus- 
tained if i?e m falls below -Re™* . We adopt a value of 
Re™ lt = 100 in this paper, following the value used in 
paper 1 and Fromang et al. (2002). We are able to cal- 
culate the distribution of Re m within our disks. Regions 
with Re m < 100 are deemed to be magnetically dead, and 
those with Re m > 100 magnetically active. 

When it comes to diffusion we consider the two ex- 
tremes with the diffusion coefficient T> — v and T> = 0. 
Models with T> = v are labelled modellD and model3D. 
Models with T> = reduce to the pure kinetic models 
modell and model3 already discussed in paper 1. 

We solved the equations by integrating over a time in- 
terval of 100,000 yrs. Hence, the ionisation fraction x[e~] 
is a function of time t, and in principle so is the location 
of the transition zone. For all kinetic models, however, the 
change in the vertical location of the transition zone at all 
cylindrical radii in the computational domain was below 
the grid resolution for t > 50, 000 yrs. 

3.1. Diffusion versus Recombination 

In order for diffusion to modify the ionisation fraction of 
material near the transition zone in protoplanetary disks, 
and hence modify the size of the dead zone, the vertical 
diffusion time should be less than the recombination time 
for free electrons. 

We consider diffusion across a scale height at R ~ 5 
AU in the heavier disk model. The scale height of this 
disk is given by H/R ~ 0.1 (see paper 1). Figure shows 
that the kinematic viscosity at 5 AU is v ~ 3 X 10 15 cm 2 
S _1 , which we adopt as the value for the vertical diffusion 
coefficient D. The diffusion time is given by 

H 2 

We consider a simple chemistry (described in detail in 
section I3.2f) involving a representative molecule and its 
ion, "m" and "m + ", and a representative heavy metal and 
its ion, "M" and "M + ". In regions well shielded from X- 
rays, the recombination rate of free electrons is given by 



625 yrs. 



(G) 



dx [e ] 
dt 



= -/j 1 x[e~]x[m+]Ar[H 2 



iM + )N[R 2 ] (7) 



where k\ is the rate coefficient for recombination between 
molecular ions and electrons, k 2 is the equivalent for heavy 
metal ions, and ^V[H 2 ] is the number density of hydrogen 
molecules. The two scenarios of interest are the metal-free 
case and the metal-dominated case. 

In the metal free case, where the electron fraction 
x[e~] = x[m + ], equation J7J becomes 



Re„ = 



He* 



(5) 



dx[e 
dt 



-fcia; 2 [e-]iV[H 2 



(8) 
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Table 1. Rate coefficients applied for the Oppenheimer & 
Dalgarno model. The values are always used as reference 
values. 



c 


Ceff 




s" 1 


a 


3 x 10" 


- 6 /VT 

-9 


cm 3 s" 1 


P 


3 x 10" 


cm 3 s _1 


7 


3 x 10" 




cm 3 s _1 



Assuming AT[H 2 ] remains constant, equation © integrates 
to give the recombination time 



1 



feiJV[H 2 ] 



1 



(9) 



where x l _ and x _ are the initial and final values of x P - 
during diffusion over a scale height. At R — 5 AU, the elec- 
tron abundance at the transition zone is a; = 5 x 10 -13 , 
as shown in figure 7 of paper 1. The rate coefficient k% = a, 
N[R 2 ] ~ 5 x 10 12 cm' 3 , and T ~ 100 K (see figures 1 and 
2 in paper 1). The recombination time is then ~ 15 
days. It is clear that diffusion is not expected to have a 
significant impact on the dead zone structure in this case. 

The rate coefficient appropriate to the metal- 
dominated case is ki = 7, which gives a recombination 
time of t~ ~ 4200 yrs. Comparing this with the diffusion 
time w = 625 yrs, it is clear that a protoplanetary disk 
with a significant population of gas-phase heavy metals 
should have the structure of the dead zone modified sub- 
stantially by the action of turbulent diffusion. 

3.2. Oppenheimer & Dalgarno model 

The underlying kinetic model of model ID involves two ele- 
ments, five species/components, and four reactions. These 
species are free electrons "e~", a representative molecule 
"m" , a heavy metal atom "M" , and their ionized coun- 
terparts "m + " and "M + ". The rate coefficients for the 
ionisation of m and for dissociative recomination of ra 
are given by £ and a, while the rate coefficients f3 and 
7 apply to the charge-transfer reaction between m + and 
M and the radiative recombination of M + , respectively. 
The rate coefficients are listed in tableland were used in 
all modellD calculations unless stated otherwise. Results 
obtained in paper 1 showed that the dead zones obtained 
using model 1 were always smaller than those obtained for 
the complex models. We find it useful, however, to anal- 
yse modell and modellD as many of the features of this 
simple model can be helpful in understanding the more 
complex model. 

The results obtained for modellD are presented in fig- 
ure |31 which shows the column density of the whole disk 
plotted as a function of radius using the solid line, and the 
column density of the active zone only using either dotted 
lines (for which T> = 0) or dashed lines (for which T> — v). 

The first thing to note is that diffusive mixing makes 
essentially no difference to the size of the active zone when 
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Fig. 3. - modell (D) - Column densities of the whole disk 
(solid line) and of the active zones (dashed and dotted 
lines) - refering to magnetic Reynolds numbers greater 
than 100 - for different values xm- While the dashed lines 
refer to simulation with T> = v, the dotted lines refer to 
simulations with T> = 0. The reference values for the rate 
coefficients are applied. The disk parameters are a = 10 -2 
and M = 10" 7 Moyr" 1 . 



heavy metals are not included in the chemistry. This is 
illustrated by the lowest line plotted in figure^ which ap- 
pears as a dash-dotted line in the figure, but is actually a 
dashed line (representing the case with diffusion switched 
on) plotted over a dotted line (representing the case when 
diffusion is switched off). The reason for this result is sim- 
ple: the recombination time when metals are absent is very 
much shorter than the turbulent mixing time. As a parcel 
of fluid moves from the active region of the disk toward 
the dead zone, free electrons will recombine with molecular 
ions rapidly, making no change to the ionisation fraction 
of the dead region. 

Increasing the heavy metal elemental abundance from 
zero to xm = 10~ 13 and switching diffusion on also makes 
no difference to the size of the dead zone. The addition of 
such a small fraction of metals does not prevent recom- 
bination between molecular ions and electrons being the 
dominant destruction process of the free electrons. This is 
illustrated by the second lowest line in figure |3 

Increasing the heavy metal fraction to xm = 10~ 12 re- 
sults in a significant change in the structure of the dead 
zone when diffusion is switched on. Interior to R < 2 AU 
there is almost no change in the size of the active zone be- 
cause recombination between electrons and molecular ions 
is still the dominant process determining the ionisation 
fraction. Beyond R > 2 AU, however, the dead zone dis- 
appears completely when diffusion is switched on, whereas 
as a dead zone remains when xm = 10 -12 and T> = 0. The 
reason for this difference is simple. The charge-transfer 
reaction between molecular ions and neutral heavy met- 
als removes most of the molecular ions beyond R > 2 
AU. This allows the recombination between metal ions 
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and electrons to become the primary means by which free 
electrons are destroyed. The recombination time for this 
reaction, however, is much longer than for recombination 
between molecular ions and electrons, and is also longer 
than the turbulent mixing time. As a parcel of fluid moves 
from the active zone toward the dead zone, the time it 
takes is shorter than the time scale over which the free 
electrons are destroyed, thus increasing the free electron 
fraction in the dead zone and allowing it to become active. 

In order to explore further the role played by the rela- 
tive diffusion and recombination time scales, and to exam- 
ine conditions under which the reaction-diffusion model 
deviates from the pure kinetic model, we have performed 
some numerical experiments in which we have artificially 
modified recombination times. We discuss these briefly be- 
low. 

Case: x u = 

In this case, model ID reduces to a three component sys- 
tem with two reactions. The molecular ion m + and the 
free electron e~ are formed and destroyed by the same 
chemical reaction. Since m + and e~ are mixed with the 
same diffusion velocity, it is simple to demonstrate how 
the mixing time scale, tt>, and the molecular ion-electron 
recombination time scale, r„, affect the ionisation fraction 
x[e~] (where Tp a 1/v and oc 1/5:). We systematically 
changed the ratio r-p/r^ by changing the rate coefficient 
a of the molecular ion - electron recombination reaction. 
We find that the recombination time scale t& has to be 
increased by three orders of magnitude before the results 
of the reaction-diffusion model model ID deviates from the 
pure kinetic model model 1. This illustrates the basic point 
that reasonable diffusion rates will be unable to change the 
size of the dead zone obtained using the Oppenheimer & 
Dalgarno model in the absence of heavy metals in the gas 
phase. 

Case: x M ^ 

We performed a similar exercise to that just described, 
but included heavy metals with an elemental concentra- 
tion of x M > 1CP 12 . This value causes the recombina- 
tion of metal ions with free electrons to dominate over 
recombination with molecular ions beyond R > 2 AU. We 
used the reference values for the rate coefficients a and (3 
listed in tabled and modified systematically the recom- 
bination rate coefficient 7 between metals and electrons 
only. Tj oc I/7 now becomes the dominant recombination 
time scale. Figure |21 shows that the addition of metals 
with x M > 10~ 12 causes a substantial difference between 
the reaction-diffusion model and the pure kinetic model. 
We find that we need to decrease the value of by two 
orders of magnitude in order that diffusion has no impact 
on the size of the dead zone. 

One can consider this issue in reverse, and ask what 
changes need to be made to in the pure kinetic model 
model 1 in order for it to give the same sized dead zone 
as model ID with T> = v and the standard value of 7. 
We find that if r=y is increased by two orders of magni- 



tude in model 1 then the resulting dead zone is the same 
as obtained in the standard modellD. For this particu- 
lar reaction network, this indicates that diffusion acts as 
an effective reduction in the dominant recombination rate 
coefficient. 

3.3. U MIST model 

The underlying kinetic model in model3D was constructed 
by extracting all species and reactions from a given species 
set of 174 species and the UMIST database containing the 
elements H, He, C, O, N, S, Si, Mg, Fe. 1965 reactions were 
extracted from the UMIST database. A full description is 
given in paper 1. Compared with the previous model of 
Oppenheimer & Dalgarno (model ID) there are now many 
pathways through which free electrons can recombinc with 
molecular ions, each with their own associated time scale. 
Without knowing a priori which ion is dominant, it is dif- 
ficult to make estimates of the relative mixing and recom- 
bination time scales. Indeed, our previous study presented 
in paper 1 suggests that there is more than one dominant 
ion that determines the ionisation fraction, complicating 
the picture further. 

An additional complication is that the introduction of 
diffusion into the kinetic model can modify the chemical 
pathways, such that comparing the products of the pure 
kinetic models with those of the reaction-diffusion models 
becomes difficult. This effect manifests itself in our study 
by changing the identities of the molecular ions which are 
dominant in the transition region between dead and active 
zones. 

We now discuss the results obtained using model3D, 
beginning with models for which the heavy metal, magne- 
sium, was was neglected from the chemical network, before 
discussing the effects that including Mg has on the ioni- 
sation fraction and resulting dead zone structure. 

Case: x Mg = 

We begin by discussing the results of model3D applied 
to the disk model defined by a = 10~ 2 and M = 
10 -7 M Q yr _1 before discussing the results obtained in the 
disk model with a = 5 x 1(T 3 and M = 1(T 8 M yr _1 . 
As a general result we find that the location of the transi- 
tion zone (which separates the active and the dead zones) 
is very similiar when comparing the results obtained for 
both disk models by applying a reaction-diffusion model 
with T> = v and T> = 0, respectively. This is illustrated 
by figure 01 in which the column density of the whole disk 
is plotted using the solid line, the column density of the 
active zone obtained using the pure kinetic model is plot- 
ted using the dotted line, and the results of the reaction- 
diffusion model are plotted using the dashed line. We see 
that model3D with T> — v produces a slightly larger active 
zone within 0.6 < R < 3 AU, but beyond this region the 
depth of the active zone is unaffected by diffusive mixing. 

We have examined the abundances of the key molec- 
ular ions in the transition region between the active and 
dead zones. Typically the ionisation balance is not deter- 
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Fig. 4. - model3 (D) - Column densities of the whole disk 
(solid line) and of the active zones (dashed and dotted 
lines) - refering to magnetic Reynolds numbers greater 
than 100 - for Xm& = 0. While the dashed line refers to 
simulations with T> = is, the dotted line refers to simula- 
tion with T> = 0. The disk parameters are a = 10 -2 and 
M = 10" 7 Moyr" 1 . 



Fig. 5. - model3(D) - Column densities of the whole disk 
(solid line) and of the active zones (dashed and dotted 
lines) - refering to magnetic Reynolds numbers greater 
than 100 - for xm s = 0. While the dashed line refers to 
simulation with T> = v, the dotted line refers to simulation 
with T> = 0. The disk parameters are a — 5 x 10~ 3 and 
M = 10~ 8 M yr- X . 



mined by a single dominant ion, but by a small number 
of ions. We find clear differences between the pure kinetic 
model and the reaction-diffusion model in those regions in 
radius where there is the greatest difference in the depth 
of the active zone between the two models. 

HCNH+ and NH^ are the most abundant ions in the 
neighbourhood of the transition zone between 1.8 < R < 
3.0 AU obtained for the pure kinetic model. It is important 
to note that the ratio between both ion concentrations 
does not change gradually in the vicinity of the transition 
zone, but instead there is a sharp transition: above the 
transition zone HCNH+ is more abundant than NH4" by 
one order of magnitude. This ratio becomes inverted just 
below the transition zone by the same order of magnitude. 
Due to the significantly shorter recombination time scale, 
electrons recombine more efficiently with NH^ than with 
HCNH+. Hence, NH+ determines the location of the tran- 
sition zone rather than HCNH+. 

By contrast, the abundance of NH^" is significantly 
lowered by applying the reaction-diffusion model with 
T> = v. Now, the more complex ion H 4 C 2 N + is two 
orders of magnitude the more abundant than the next 
most abundant ion CH3OHJ. In addition, since mixing 
tends to reduce the concentration gradients of each species 
with height, z, the ratio between both ion concentrations 
a;[H 4 C 2 N + ]/a;[CH30H^] changes gradually when passing 
through the transition zone. In the same region between 
1.8 < R < 3 AU we found significant differences in the 
distribution of most of the species obtained for a reaction- 
diffusion model and the corresponding pure kinetic model. 
This is a clear indication that mixing may change the ki- 
netics significantly. 

We now consider the a = 5x 10~ 3 , M = 10~ 8 M Q yr _1 



disk model. The results for this model are shown in fig- 
ure We first note that this disk model generates an 
mtrmisically deeper dead zone in the outer regions of the 
disk beyond R > 3 AU, even though the ionisation rate 
due to X-rays is higher due to the smaller surface density. 
This is a temperature effect that arises because the viscous 
heating is reduced because of the smaller a and M values. 
The recombination rates scale according to k cx 1/y/T, 
leading to a lower ionisation fraction. 

We focus in particular on locations 0.6 < R < 1.1 AU 
where the most significant differences in the vertical depth 
of the active zone appear when comparing the reaction- 
diffusion model and the corresponding pure kinetic model. 

The most abundant ions produced by the pure ki- 
netic model in the neighbourhood of the transition zone 
are H 4 C 2 N + and NH 4 . These species are similarly abun- 
dant, but the recombination time for electrons with NH 4 
is much shorter than with H 4 C 2 N + , so NH4" determines 
the position of the transition zone. The corresponding 
reaction-diffusion model again results in £[NH 4 ] drop- 
ping significantly so that it no longer controls the ionisa- 
tion fraction in the neighbourhood of the transition zone. 
Instead, H 4 C 2 N + is the most abundant molecular ion, be- 
ing one order of magnitude more abundant than the next 
most abundant ion H 3 + . This is again evidence that the 
presence of turbulent mixing in a protoplanetary disk can 
change the reaction kinetics. This leads to the intriguing 
possibility that such changes may in the future act as ob- 
servational diagnostics for the presence of turbulence in 
disks. 

Case: x Mg / 

We now consider the effect that including magnesium has 
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Fig. 6. - model3 (D) - Column densities of the whole disk 
(solid line) and of the active zones (dashed and dotted 
lines) - refering to magnetic Reynolds numbers greater 
than 100 - for different values of a?Mg- While the dashed 
lines refer to simulations with T> = u, the dotted lines 
refer to simulations with T> = 0. Especially, for a reaction- 
diffusion model with 2? c ff = v no dead zones are observed 
above values x Mg > 10~ 8 . The disk parameters are a = 
10~ 2 and M = 1O- 7 M yr" 1 . 

on the models. We begin by discussing the results obtained 
for the disk model with a = 10~ 2 , M = 10~ 7 Moyr" 1 . 
The results for model3D are shown in figureHjl The dotted 
lines show the column density of the active zones obtained 
by the pure kinetic model for various values of x-^g ■ The 
dashed lines show the column density of the active zones 
obtained by the reaction-diffusion model. In general we 
find that mixing affects the size of the active zone sig- 
nificantly as soon as metals with concentrations above a 
threshold value are involved. 

The addition of magnesium with elemental concentra- 
tion x Mg = 10~ 12 makes very little difference to the size of 
the dead zone, whether diffusion is switched on or not, as 
seen by comparing figures 0] and [|J] This is because a low 
value of x Mg does not allow Mg + to become the dominant 
ion via charge-transfer reactions. 

Increasing the magnesium elemental abundance to 
^Mg — 10~ 10 changes the picture dramatically when dif- 
fusion is switched on. The pure kinetic models with these 
magnesium abundances included contain no regions in the 
disk beyond R > 0.4 AU where the disk is fully active, al- 
though the dead zone tends to shrink beyond R > 2 AU. In 
this region Mg + followed by HCNH+ are the most domi- 
nant ions, while in the region 0.4 < R < 2 AU the NH|" ion 
is dominant. The reaction-diffusion model model3D with 
T> = v and a; Mg = 10 -10 produces a dead zone of similar 
depth to the pure kinetic model between radii 0.4 < R < 2 
AU. The identity of the dominant molecular ion in this re- 
gion changes, however, when diffusion is switched on. The 
abundances of HCNH+ and NH4" are found to decrease 
significantly, and are replaced by the ion H4C2N+. The 



Fig. 7. - model3(D) - Column densities of the whole disk 
(solid line) and of the active zones (dashed and dotted 
lines) - refering to magnetic Reynolds numbers greater 
than 100 - for different values of iMg- While the dashed 
lines refer to simulations with T> = v, the dotted lines 
refer to simulations with T> — 0. The disk parameters are 
a = 5 x 10" 3 and M = 10~ 8 Moyr" 1 . 

disk is rendered fully active beyond R > 2 AU where Mg + 
is the dominant ion and the recombination time increases 
above the mixing time. 

model3D with T> = v and x M<y = 10~ 8 produces a fully 
active disk, in contrast to the corresponding pure kinetic 
model model3. For these higher magnesium abundances, 
Mg + becomes the dominant ion through the charge- 
transfer reaction with molecular ions, and the abundance 
of NH^" is much reduced. The long electron recombina- 
tion time associated with Mg + allows diffusive mixing to 
maintain a disk that is fully active. 

We now consider the disk model with a = 5 x 10~ 3 , 
M — 10~ 8 M©yr _1 . The results are presented in figure[7| 
where the column density of the active zones obtained us- 
ing the pure kinetic models are plotted using dotted lines. 
Dashed lines represent the column density of active zones 
obtained using the reaction-diffusion model. 

Adding magnesium with an elemental abundance 
£Mg = 10~ 12 results in very little change in the size of the 
dead zone whether diffusion is switched on or not (com- 
pare figures and EJ). 

Increasing the magnesium elemental concentration to 
XMg = 10~ 10 causes substantial changes to occur when 
diffusion is switched on. Interior to R < 0.7 AU there is 
no change because the molecular ions NH4" and H4C2N+ 
are dominant over Mg + there, both with and without dif- 
fusive mixing. Between 0.7 < R < 3 AU the disk becomes 
fully active when diffusion is switched on, but maintains 
a significant dead zone there in the pure kinetic model. 
The influence of Mg + allows the recombination time to 
be increased so that mixing can render the disk active. 
Beyond R > 5 AU we obtain the surprising result that 
diffusive mixing causes the depth of the dead zone to in- 
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crease rather than decrease. The reason for this is that the 
ionisation fraction at large radii in the pure kinetic model 
is always close to the critical value for the disk to be ac- 
tive. The addition of diffusion appears to have a diluting 
effect on the electron fraction in the active zone, which 
causes the dead zone to increase in size. 

Increasing the magnesium elemental concentration to 
^Mg = 10~ 8 causes further changes to occur when diffu- 
sion is switched on. The Mg + ion is dominant throughout 
the disk, and the disk model is rendered entirely active 
between the radii 0.4 < R < 6 AU. We again find that 
the dead zone becomes deeper at radii beyond R > 8 AU 
when diffusive mixing is effective. The reason is the same 
as given above: the ionisation fraction at these radii is al- 
ways close to the critical value required to render the disk 
active. Diffusive mixing dilutes the electron fraction there, 
producing a deeper dead zone. 



4. Summary 

We have presented the results of calculations that exam- 
ine the ionisation fraction in protoplanetary disks mod- 
els. Our models compute the chemical evolution of the 
gas, with X-rays from the central star being the primary 
source of ionisation, and also include vertical diffusion of 
the chemical species. This diffusion is designed to mimic 
the effects of turbulent mixing arising from the MHD tur- 
bulence that is thought to drive accretion in protostellar 
disks (e.g. Balbus & Hawley 1991). 
The main findings of this work are: 

— The simple Oppcnhcimer & Dalgarno (1974) kinetic 
model produces smaller dead zones than a more com- 
plete kinetic model based on the UMIST database. 
This is true whether diffusionis switched on or not. 

— All models which did not include heavy metals (magne- 
sium) give rise to substantial dead zones. The inclusion 
of diffusion in these cases makes very little difference 
to the size and depth of the dead zone. This is because 
the recombination of molecular ions with electrons oc- 
curs too rapidly for turbulent mixing to be effective. 

— The addition of heavy metals (magnesium) to the re- 
action networks can give rise to dramatic changes in 
the sizes of dead zones when diffusion is included. This 
is because the recombination time between metal ions 
and electrons is orders of magnitudes longer than that 
between molecular ions and electrons, and is longer 
than the diffusive mixing time. The addition of met- 
als, combined with charge-transfer reactions, removes 
most of the molecular ions and replaces them with the 
longer lived metal ions. 

— For the Oppcnhcimer & Dalgarno (1974) network, the 
addition of heavy metals with fractional abundance 
xm = 10~ 12 leads to the dead zone disappearing be- 
yond R > 2 AU when diffusion is switched on. In the 
absence of diffusion a substantial dead zone exists in 
this region of this disk. 



— For the model based on the UMIST database, the ad- 
dition of magnesium to the chemical network with a 
fractional abundance of XMg = 10~ 10 leads to the dead 
zone disappearing beyond R > 2 AU when diffusion is 
switched on. When the abundance of magnesium is in- 
creased to XMg = 10 the dead zone disappears com- 
pletely. Models in which diffusion is not present, but 
with the same magnesium abundances have substan- 
tial dead zones. 

— We find that the changes in the size of dead zones 
when diffusion is switched on do not only arise be- 
cause free electrons are mixed into the dead zone, but 
also because the chemical pathways are modified by 
diffusion. One effect of this in our models is that the 
dominant ions near the transition zone between dead 
and active zones change their identity when diffusion 
is initiated. This raises the interesting possibility that 
there may be chemical signatures of turbulent diffusion 
that can act potentially as observational diagnostics of 
turbulence in protoplanetary disks. 

— For some disk models, the inclusion of diffusion can 
cause the dead zone depth to increase in the outer 
regions of the disk. This arises because these par- 
ticlar models predict that the electron fraction is only 
slightly above the critical value for the disk to be ac- 
tive in these regions when diffusion is switched off. 
Switching diffusion on can dilute the electron fraction 
in these regions, leading to a deeper dead zone. 

Overall, our models indicate that turbulence in protoplan- 
etary disks can be a self-sustaining process, providing that 
certain criteria are met, and indicate that dead zones can 
be reduced or removed altogther through turbulent mix- 
ing processes. These criteria are: 

(z) There are sufficient magnesium atoms available in the 

gas phase so that recombination between magnesium ions 

and electrons becomes the dominant process by which the 

local ionisation fraction is determined. 

(ii) The turbulent diffusion time scale is shorter than the 

dominant recombination time on which free electrons are 

removed. 

A potentially serious omission from our models is the 
existence of small dust grains. These are known to sweep 
up free electrons (and metal atoms) very efficiently, and 
when included with the interstellar abundance they sub- 
stantially increase the size and depth of dead zones in 
protoplanetary disks (see Sano et al 2000; paper 1). It is 
very likely that these grains will need to be substantially 
depleted by grain growth before the effects of turbulent 
mixing described in this paper are realised, as the removal 
rate of electrons from the gas phase is too large to allow 
diffusion to be effective. 

An interesting question is whether a disk that has a 
substantial dead zone, and which is subject to a tempo- 
rary period of intense ionisation (by large X-ray flares of 
the type recently report by Feigclson et al. 2005, for ex- 
ample), can enter into a self-sustaining turbulent state 
in which the turbulent mixing maintains a fully active 
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disk. This question is currently under investigation and 
will be the subject of a future publication. An equally 
interesting question is whether a disk with a substantial 
dead zone, sandwiched between active zones near the disk 
surface, can become globally active due to turbulent mo- 
tions overshooting into the dead zone and transporting 
free electrons into it. This question cannot be addressed 
by calculations of the type presented in this paper, but 
their results suggest that this may be possible under cer- 
tain circumstances. 
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Appendix A: 

The diffusion flow Jf of the i component may be ex- 
pressed in various ways depending which reference veloc- 
ity v a is used. The molar diffusion flow Jf is 



Jt = n t (v, - v Q ) . 



(A.l) 



The reference velocities v a are defined as weighted aver- 
ages of the velocity v ; of the i th component 



v a = ^a 4 v 4 , |Vfl, = 1 



(A.2) 



with a i as normalised weights. 

Since the mass is conserved, the diffusion fluxes of the s 
components are not independent. Specifying the reference 
velocity v a by the mean molar velocity v m 



E 

i=l 



(A.3) 



where x i denotes the molar fraction of the i component 



(A.4) 




the molar diffusion fluxes are constrained by 



5> 



0. 



(A.5) 



There are several contributions to the diffusion flux. We 
assume that the most important contribution to the diffu- 
sion flux is caused by concentration gradients rather than 
due to pressure gradients, external forces, and coupled 
effects (cross-effects). Diffusion caused by concentration 
gradients is called "ordinary" diffusion. Considering or- 
dinary diffusion processes only, we often use the phrase 
"diffusion" instead of "ordinary" diffusion. The (ordinary) 
diffusion is described by the phenomenological equation, 
which is well-known as Fick's first law for a binary system 
s = 2. By analogy to the molar diffusion flow for a binary 
system, we define an effective molar diffusion flow JJJ i 
for the diffusion of the i component in a mixture by 



(A.6) 



The effective binary diffusivity T> cS i for species i in a mix- 
ture can be calculated from the diffusivity T>^ of the pair 
i—j in a binary mixture by applying the so-called "Stefan- 
Maxwell equation" . Especially, for systems in which all the 
T>^ are the same T> eS i is now given by 



V 



eS,i 



(A.7) 



Because of the definition of the molar diffusion flow Jf , see 
equation l[A.ljl . the velocity of the i component and the 
reference velocity v a are related to eachother. In particu- 
lar, in a molar diluted system 1 we deal with, one can ap- 
proximate the reference velocity by v" = v m = 0. Similiar 

1 i.e. for m ; <C 1, £ = 1, 2, . . . , s — 1, n s ~ 1 
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applies for mass diluted systems with v" = v n = where 
v„ denotes the n th component velocity. This is in agree- 
ment with the hydrostatic equilibrium assumption in z 
direction. Finally, equation @ can be simplified by: 

~at = ~d~z ( nV ^d^ Xi ) + 51 v u J j> (* = ■ ■ ■ ' s ) ( A - 8 ) 

constraint by equation (|A.5(1 . Note that we simplified the 
symbol for the effective binary diffusivity for species i in 
a mixture by using T> instead of T> cS i . 

Because of the molar diluted system we consider and 
due to the fact that the most abundant components 
(which are molecular hydrogen and helium) do not con- 
tribute to the molar diffusion fluxes, the mass is conserved 
even by neglecting equation (|A.5I) . 



